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We present results on spatio-temporal correlations in the so-called mean drag version of the Durian 
bubble model in the limit of small, but finite, shearing rates, 7 . We study the rheology, diffusion, 
and spatial correlations of the instantaneous velocity field. The quasi-static (QS) effective diffusion 
co-efficient, D e , shows an anomalous system size dependence indicative of organization of plastic 
slip into lines along the directions of maximum shearing. At higher rates, D e decays like 7 -1 / 3 . 
The instantaneous velocity fields have a spatial structure which is consistent with a set of spatially 
uncorrelated Eshelby transformations. The correlations are cut off beyond a length, f. £ ~ '-y -1 / 3 
which explains the D e ~ j- 1 / 3 behavior. The shear stress, a, follows a similar rate dependence with 
5a = a — a y ~ "y 1//3 where a y is the yield stress observed in the QS regime. These results indicate 
that the form for the viscous dissipation can have a profound impact on the rheology, diffusion and 
spatial correlations in sheared soft glassy systems. 


In recent years, it has become clear that shear flow 
in many types of amorphous solids - structural glasses 
such as metallic or polymer glasses; soft glasses such 
as pastes, emulsions, or foams; sheared granular mat¬ 
ter - is governed by localized plastic shear transfor¬ 
mations W- These solids behave like elastic bod¬ 
ies below a yield stress, a y , but flow like viscous liq¬ 
uids above. At low shearing rates, the response be¬ 
comes bursty and intermittent, resembling other slowly 
driven out-of-equilibrium systems such as pinned elas¬ 
tic manifolds, Barkhausen noise in disordered ferromag- 
nets, martensitic phase transformations, or dislocation- 
mediated crystal plasticity |jy. One central question is 
how finite driving rate affects the spatio-temporal fluctu¬ 
ations and how this impacts the average response. 

In the present case of amorphous solids, models have 
focused at the particle scale [2] |3] I6HT6] and at coarser, 
meso-scale levels nzmsnt At the meso-scale, one en¬ 
visions plastic flow occurring due to localized yield¬ 
ing events that are observed in experiments and in 
the particle-scale models. One may construct mean 
field theories based on these localized yielding events. 
Some examples are the shear transformation zone (STZ) 
theory pi;, the soft glassy rheology (SGR) mo del [32] [33], 
and the Hebraud-Lecquex model (HL) 31i 35 . These 
models vary in their assumptions about the distribution 
of states, the barriers seen by the states, and the dynam¬ 
ics of the yielding events, but they all predict a so called 
Herschel-Bulkley (HB) power-law rate dependence of the 
shear stress where 5a ^ 7^. However, the exponent, /?, 
varies from model to model. 

One can, alternatively, construct explicit real-space 
models. These so-called elasto-plastic models (EPMs) 
were pioneered by Bulatov and Argon and many varia¬ 
tions have been proposed in the intervening years [nr 
EDI. In most, with a few exceptions mm, the loads 
are transferred instantaneously across all space as the 
site in question undergoes a yielding event. In one of the 
few EPMs where loads are transferred dynamically, a 
5a ~ 7 0 5 rheology was also observed [20]. Recently, Liu 


et. al. [36] . have shown that in a rate dependent EPM, 
the rheology exhibits a crossover from a non-trivial uni¬ 
versal scaling regime at low rate where 5a ~ 7 0 65 to a 
mean field behavior at higher rate where 5a ~ 7 0 51 . 

Also at the particle scale, one generally observes HB 
rheology, often with a 1/2 exponent. Molecular dynam¬ 
ics (MD) simulations on low temperature Lennard-Jones 
(LJ) glasses [12] are in agreement with the HL model 
which predicts 5a ~ 7 1 / 2 . In these simulations, there 
was also a power-law dependence on rate in the corre¬ 
lation length governing the plastic strain field and the 
diffusion coefficient, and a connection between these spa¬ 
tial correlations and the HB exponent was suggested. In 
foam simulations using Durian’s bubble model with a lin¬ 
ear drag imposed on relative velocity of particles m and 
more sophisticated models incorporating non-linear elas¬ 
ticity and non-linear drag between the particles [38] along 
with corresponding experiments on similar systems [38L 
[40] all show behavior consistent with 5a ~ 7 1 / 2 . 

Between the particle-based simulations [12] [37] [38] , the 
EPMs [20] [36] , and the HL model [34] [35] , one might ex¬ 
pect 5a ~ 7 1 / 2 to be universally applicable to many sys¬ 
tems where the basic picture of local yielding and stress 
redistribution holds. We show here that a simple variant 
of Durian’s bubble model mm falls outside this class 
with an HB exponent of roughly 1/3. Furthermore, the 
class of experiments which are closest to this model - 
amorphous Bragg-Nye bubble rafts - mm also show 
a HB rheology with an exponent consistent with 1/3. 

In this work, we study a version of the bubble 
model m m where the drag on the particles arises 
from motion with respect to the background suspending 
fluid; the so-called mean field drag (MFD) variant. A 
mean-drag term could be appropriate for modeling sys¬ 
tems such as particles at an air-water interface [43H47] 
where dissipative forces on a particle may be governed by 
the generated subsurface flow. Early results on the MFD 
model were restricted to small systems and large shear¬ 
ing rates mm- More recent studies have focused on 
the behavior near the jamming transition. Here, we work 
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with much larger systems and at much lower rates than 
previous studies and at volume fractions well above the 
random close packing point (precisely the same volume 
fractions and particle size and stoichiometry as studied 
by Ono et. al. [41J [42] ). 

We show that the MFD bubble model gives an HB ex¬ 
ponent of roughly 1/3. Surprisingly, there are long range 
correlations present in the instantaneous velocity fields 
that are of the classical Eshelby form one would expect 
for the displacement fields induced by a local shear trans¬ 
formation in an elastic matrix [48]. The Eshelby corre¬ 
lations are cutoff beyond a length, £, which scales like 
^ ^ 7 -1 / 3 . £ governs the effective diffusion co-efficient, 
D e , as in the MD simulations [12], with D e ~ £, and 
the HB exponent via 8a ~ l/£. This gives an effective 
Stokes-Einstein relation: D e 8a = const as in the MD 
case [ 12 ] (where we define 5rj = 8a j 7 ). 

We consider two dimensional ( 2 D) system of soft 
disc particles in a bi-disperse mixture of equal num¬ 
ber of small and large discs with diameter ratio, D$ : 
D l = 1 : 1.4 j4TJ 02]. We simulate this mixture at 
zero temperature using the bubble model introduced by 
Durian|49j 50,. The harmonic interaction between discs 
i and j is, Uij = k e 8 ^/2 if 8{j < 0 and zero otherwise, 
where Sij = 2 77 j /(Di + Dj) — 1 , is the overlap distance, k e 
is the elastic spring constant. Viscous dissipation is taken 
into account in a mean-field fashion with the total drag 
force on disc i, FP = —b(vi —yi 7 #), where b is the damp¬ 
ing parameter, 77 is the location of the particle projected 
along the flow-gradient direction, x is the unit vector in 
the flow direction, and 7 is the imposed shearing rate. 
Since the particles are considered massless, the dynamics 
are given by balancing elastic and drag forces. The only 
timescale in the model is the timescale for elastic relax¬ 
ation under the drag force: td = f; in our simulations we 
have used 6=1 and k = 1. We report all lengths in units 
of Ds and all times in tr. We define the volume frac¬ 
tion, </>, as ^(TVlD 2 + NsD 2 s )/AL 2 , where Nl and Nr are 
the number of big and small particles respectively, and 
L is the simulation box size. We set cj> = 0.9 which is far 
above the jamming point, <fij ~ 0.843 in two dimension. 
We report results for different sizes, L = 40,80,160 cor¬ 
responding to a total number of, N = 1240,4960,19840 
particles. Lees-Edwards boundary conditions are used 
to impose simple shear flow [51] . We refer to the flow 
direction as x and the gradient direction as y. We use 
the standard Irving-Kirkwood expression for the shear 
stress a = (1 /L 2 ) JW F x ijr y ij [51] [52]. In this letter, we 
report on the steady state achieved beyond 50% strain, 
for different L and rates, 7 , and all statistics are taken 
between 50% and 150% strain. All velocities discussed 
below are the non-afRne velocities defined with respect 
to the background flow. 

In Fig. Ef> we plot the effective diffusion constant D e = 
hmA 7 -s>oo(A 2/ 2 )/267 vs. 7 for the three different system 
sizes, L = 40, 80 and 160. Ay is the particle displacement 




FIG. 1: (a) D e vs. 7 for systems of size L = 40, 80, and 160; 
thin dashed line has a slope of —1/3. Inset: D e /L vs. Ly 1 ^ 3 . 

(b) a vs 7 ; thin dashed line has a slope of 1/3. Inset: effective 
athermal Stokes-Einstein temperature, D e 8a , for the L — 160 
system for a y = 0 . 0011 . 

in the direction transverse to the flow for a time during 
which a strain of amplitude, A 7 = 7 At, was applied. At 
large rates, D e , is independent of L and follows a 7 -1 / 3 
power law with remarkable precision and for over four 
decades of rate for the case of the L = 160 system. At the 
lowest rate in the QS regime, a plateau is clearly visible 
for L = 40 and 80. In the inset, we plot D e /L against 
L 7 1 / 3 . The data shows a good collapse, indicating that 
the QS diffusion, Dqs ~ L. This behavior has been 
understood to arise from system spanning lines of slip 
in MD simulations of LJ glasses [7j 111413] . The bubble 
model exhibits these same slip lines in the QS regime. 
It is, rather, the rate dependence which is different here 
- both in the rheology and diffusion. Because of the 
relatively weak rate sensitivity, the L = 160 system shows 
a strongly rate dependent D e even at 7 = 10 -7 , as it has 
not yet reached its QS plateau. 

In Fig. Et, we plot cr vs. 7 . At low rates, one ap¬ 
proaches the yield stress, a y . At higher rates, one ap¬ 
proaches a power-law regime which is described reason¬ 
ably well by 8a ~ 7 1 / 3 for over a decade. We cannot 
rule out a cross-over to a different behavior at very low 
rates as one approaches the QS limit. Regardless, the 
rate sensitivity is significantly less than 7 1 / 2 . Ono et. 
al. [42] have already studied the diffusion explicitly and 
the rheology implicitly (via the mean-squared-velocity); 
but here we explicitly demonstrate the 7 ±:L / 3 scaling and 
the finite size effects in the QS regime. In the inset, we 
plot D e 8a for the L = 160 system. D e 8a can be consid¬ 
ered a kind of effective Stokes-Einstein temperature. Our 
results for D e a agree with Ono et. al. for the 7 where 
the studies overlap, but we show here that one can ob¬ 
tain a remarkably constant value of D e 8a for over three 
decades in rate for the L = 160 system. 

In Fig. [2] we plot a typical snapshot of the y- 
component of the particle velocities at two typical rates, 
(a) 7 = 10 -3 and (b) 7 = 10 -5 , and their respective 
time-averaged spatial autocorrelation functions, C Vy , in 

(c) and (d). As usual, we define C Vy (R ) = (v y (R + 
r,t)v y (r,t))(ftjy The large coherent spatial structures 
are obvious. These correspond to long lines of particles 
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FIG. 2 : Top: A typical map of l^v y for (a) 7 = 10 3 
and (b) 7 = 10 -5 using L — 160. Bottom: Autocorrelation 
of y -velocities C Vy (R)/C Vy (x = 1) for (c) 7 = 10 -3 and (d) 
7 - 10 - 5 . 


L = 160 system is just starting to show system-size de¬ 
pendent behavior at the lowest rate. At higher rates, the 
data is well described by a £ ^ 7 -1 / 3 power law. There 
are deviations from scaling in the high 7 , small £ regime, 
below about £ ~ 5. Nonetheless, we can observe over a 
decade of scaling for the L = 160 system. In the inset, 
we plot £/L vs L 7 1 / 3 showing the £qs ~ L quasi-static 
scaling [53] . 



FIG. 4: (a) A typical map of 10 4 e for 7 = 10 5 ,L = 160. 

(b) Ce(R — x)/Ce(x — 2 ) for different 7 , L = 160; thin 
dashed line has a slope of - 2 . Inset: 10 2 Ce(R)/Ce(r = 0) for 
7 = 10 -5 . 


which are all slipping together at the same instant. The 
sharp jumps from blue to red as one traverses the image 
from left to right indicate discontinuities in the veloc¬ 
ity field with counter-clockwise vorticity. The real-space 
correlation functions reflect the visual impression. The 
//-velocities have strong correlations along the y direc¬ 
tion. They are correlated along the x-direction for some 
distance and eventually become anti-correlated. One 
would naturally associate the cross-over from correlation 
to anti-correlation with the typical spacing between the 
slip lines. 




FIG. 3: (a) C Vy (R — x)/C Vy {x — 1 ) for different 7 , L = 

160; black dashed line shows zero correlation, (b) £ vs 7 for 
different L. Inset: £/L vs L 7 1 / 3 ; thin dashed line has a slope 
of — 1 . 

In Fig. IT we plot the traces of C Vy , along the x- 
separations, normalized by the x = 1 values, for various 
shearing rate for L = 160. In Fig. [ 3 J 3 , we plot the loca¬ 
tion, £, of the minima of each C Vy curve as a function of 
rate for various system size. As with the D e plots, for 
L = 40 and 80, we see a clear QS plateau at the lowest 
rates where £ saturates near the system size, while the 


In Fig. IT we plot a typical snapshot of the sym¬ 
metrized shear-strain-rate field, e A ( d x v y + d y v x )/2 , for 
7 = 10 -5 for precisely the same configuration as shown 
in Fig. |2j>. [54] The vertically oriented slip lines with 
counter-clockwise vorticity, visible as sharp blue-red dis¬ 
continuities in Fig. [2 }d, appear now as red vertical lines. 
Horizontally oriented slip lines with clockwise vorticity, 
which would have been impossible to detect in the v y 
fields, appear here as well. In real space, these slip lines 
appear to form a grid with a characteristic spacing con¬ 
sistent with £. 

In Fig. [ 4 ) 3 , in the inset, we plot the respective time- 
averaged spatial autocorrelation function. The angular 
variation has a quadrupolar symmetry reminiscent of 
the strain fields one would obtain from Eshelby inclu¬ 
sions g8] with strong correlations along the directions of 
maximum shear and strong anitcorrelations 45 degrees 
away. We have checked that the angular dependence, 
at intermediate distances, is consistent with the cos A0 
dependence one would observe for an Eshelby transfor¬ 
mation. 

In Fig. [ 4 ) 3 , in the main plot, we plot the correlations 
along the directions of maximal shear, arbitrarily nor¬ 
malized to unity at x = 2, for L = 160. For the highest 
rates with £ < 5, corresponding to the regime with de¬ 
viations from scaling observed in Fig. [3] the curves have 
a non-universal form. At lower rates, where particle-size 
effects are less important, one observes a C ~ r~ 2 power 
law for x < £. The r ~ 2 power law is precisely what one 
would obtain for uncorrelated Eshelby transformations. 

In summary, we have shown that the bubble model, 
with so-called mean field drag, exhibits surprising corre¬ 
lations in the instantaneous velocity fields. These cor- 













































4 


relations are of precisely the Eshelby form in the limit 
of low shearing rate. Strong correlations - if not pre¬ 
cisely of the Eshelby form - were to be expected at low 
shear rate in the finite time displacements, in analogy 
with MD simulations, but it was surprising to see them 
here in the velocity field and surprising to see them in 
precisely the Eshelby form. Although we find essentially 
the same connection between the correlation length and 
effective diffusion co-efficient as the MD simulations of 
Lemaitre and Caroli [ 12 ], both the rate dependence of 
the correlation length and the connection between that 
length and the rheology is qualitatively different here. 

In [ 12 ], a length, l, was introduced to rationalize 
the rate dependence of D e . The relation D e ~ £ (for 
£ << L), remains the same here and can be understood 
as a simple kinematic consequence of plastic deformation 
essentially organizing along lines as in ref [12]. However, 
the rate dependence of £ is different: here it scales like 
7 -1 / 3 whereas, in reference [12], l ~ 7 -1 / 2 . The argu¬ 
ment in [ 12 ] for £ ~ 7 -1 / 2 was quite general and relied 
essentially only on: i) the assumption of a 7 indepen¬ 
dent timescale, r, for elementary shear transformations 
and ii) a basic picture of elementary shear transforma¬ 
tions releasing a characteristic, 7 independent, strain (or, 
equivalently, stress). It will be important in the future 
to try to understand precisely why this argument breaks 
down for the MFD bubble model. 

Furthermore, the connection between £ (or, equiva¬ 
lently D e ) and 8 a is different than in m- Lemaitre and 
Caroli argued that a characteristic time for stress relax¬ 
ation, r, scaling linearly with the correlation length along 
with the observed £ ~ 7 -1 / 2 behavior could explain the 
8a ~ 7 1 / 2 rheology. Here, since we observe £ ^ 7 -1 / 3 
and 8a ~ 7 1 / 3 , we would need to invoke a r ~ £ 2 rather 
than r ~ £ x relation to explain the connection . 

Interestingly, we observe a rate independent effective 
athermal Stokes-Einstein relation connecting the diffu¬ 
sion to the rheology: D e 8a = D8r] = const. Note that 
our results agree with Ono et. al. m for D e a, but 
we find it interesting that D e 8cr remains constant over 
a much broader range of shearing rates than D e a. Of 
course, at larger rates, 7 > 10 -2 , where £ « 0 ( 1 ) one 
would expect an increase in D e a , as observed in [42], but 
we are not interested in this regime for purposes of this 
study. 

It seems to us that the emergence of a rate indepen¬ 
dent effective Stokes-Einstein temperature is a more com¬ 
pelling - and apparently general - connection between 
the diffusion and rheology than a particular relationship 
between correlation length and stress relaxation time as 
in [ 12 ]. We further speculate that D e 8a may be related 
to various effective temperatures which arise in the STZ, 
SGR, and HL models. In particular, Langer argues eh. 
on very general grounds, that the effective temperature 
should become rate independent at low rate whenever 
the timescale of the external drive is much longer than 


any internal relaxation times. This rate independent T e ff 
would then give the connection between the diffusion and 
rheology via an effective athermal Stokes-Einstein rela¬ 
tion mm- 

In addition to the questions about the origin of the 
rate dependence of £ and the relation between £ and 8a, 
there are some other obvious open questions, i) If one 
imposes a drag on local shearing rates - as in the orig¬ 
inal version of the bubble model - instead of velocities, 
does one observe the same rate dependence on 7 of 8a, 
D e , and £? We are in the process of conducting these 
studies now. ii) Does one recover the same scaling laws 
when physical stress propagation - in terms of the mag¬ 
nitude and type of drag - is implemented in the EPMs? 
We note that Jagla’s model m , which includes physical 
stress propagation, in its current form, damps the local 
strain rate, as in the original formulation of Durian’s bub¬ 
ble model, rather than the non-affine velocity, as in the 
present MFD version. It would be quite interesting and 
relatively straightforward to modify this model to imple¬ 
ment velocity damping rather than strain-rate damping 
to see if one would obtain the y ^ 1 / 3 scaling laws, iii) Do 
the y ^ 1 / 3 scaling laws hold in 3D? These are some of the 
questions which will need to be addressed on the road 
to understanding the rate dependent response of various 
athermal, sheared amorphous materials. 
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